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The low-level interpretation of images provides constraints on 3D surface shape at multiple reso¬ 
lutions, but typically only at scattered locations over the visual field. Subsequent visual processing 
can be facilitated substantially if the scattered shape constraints are immediately transformed into 
Risible-surface representations that unambiguously specify surface shape at every image point. The 
required transformation is shown to lead tej an ill-posed surface reconstruction problem. A well- 
posed variational principle formulation is obtained by invoking “controlled continuity,'’ a physically 
fionrestrictive (generic) assumption about surfaces which is nonetheless strong enough to guarantee 
Unique solutions, flic variational principle, which admits an appealing physical interpretation, is 
Ideally discretized by applying die finite element method to a piecewise, finite clement represen¬ 
tation of surfaces, lliis forms the mathematical basis of a unified and general framework for 
computing visible-surface representations. The computational framework unifies formal solutions 
to the key problems of (i) integrating multiscale constraints on surface dcpdi and orientation from 
njiultiplc visual sources, (ii) interpolating diesc scattered constraints into dense, piecewise smooth 
surfaces, (iii) discovering surface depth and orientation discontinuities and allowing them to restrict 
interpolation appropriately, and (iv) overcoming the immense computational burden of fine resolu¬ 
tion surface reconstruction. An efficient surface reconstruction algorithm is developed. It exploits 
multi resolution hierarchies of cooperative relaxation processes and is suitable for implementation on 
massively parallel networks of simple, locally interconnected processors. The algorithm is evaluated 
empirically in a diversity of applications. 
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l. Introduction 


Over thirty years ago, J.J. Gibson [1950] made the seminal conjecture that natural human perception 
amounts to the perception of visible surfaces. 'Hie explicit representation of visible surfaces, an 
intermediate goal of computational vision, has since attracted considerable interest. 

The computational framework offered in this paper addresses, in a unified way, certain visual 
information processing tasks involved in the representation of visible surfaces. Particular emphasis 
is placed on utilizing highly parallel, cooperative processing to integrate surface shape information 
over multiple visual sources, to fuse it across a multiplicity of spatial resolutions, and to maintain the 
global consistency of the resulting distributed shape representations. The issues are first investigated 
in terms of a surface reconstruction model rooted in mathematical physics. This formal analysis is 
augmented by an empirical study of the resulting algorithms, which feature multiresolution iterative 
processing within hierarchical surface shape representations. The approach is guided by current 
knowledge of how humans perceive visibly surfaces, while applications in machine vision provide 
a testbed for the algorithms. 

The remainder of this introductory section examines the role of surface representations in early 
visual processing, outlines the key computational problems that will be of primary concern, and 
reviews some relevant prior work. 

1.1. Early Visual Processing and Visible-Surface Representations 

Early vision comprises a set of processes which specialize in recovering the physical properties of 
visible surfaces in a 3D scene from 2D images of the scene. They apply generic assumptions about 
the physical world and the imaging process to infer 3D surface shape constraints by interpreting 
specific image cues, such as stereoscopic disparity, motion, texture, contours, and shading. These 
conceptually independent shape estimatiorj processes fall into two broad categories. 

The first category comprises what arc commonly referred to as correspondence processes. They 
operate over multiple image frames of a scene taken across space or over time. Paradigm examples 
are stereopsis and structure from motion (see, e.g., the review articles [Poggio and Poggio, 1984] 
and [Uliman, 1983]). Stereopsis is driven by computations on typically two image frames taken 
simultaneously, but from different spatial positions, flic basic structure from motion computation 
involves frames taken from the same position, but at different times. If correspondences can be 
established across the frames, between image features which originate from the same point on 
a visible surface (not a trivial problem), then the depth (i.e., 3D distance) to such points can 
be estimated by triangulation, given the disparity (i.e., 2D displacement) between corresponding 
features as well as some knowledge of the imaging geometry. 

The second category of shape estimation processes involve computations on a single static 
frame. Perspective projection of 3D scenes onto images imparts a systematic distortion to imaged 
surface properties such as shading, texture, and contours. A major part of this distortion can be 
attributed to the relative orientations of visible surfaces with respect to the viewer. In principle, 
it is possible to estimate surface orientation by measuring and interpreting such distortions in the 
image. This is the basis of practical approaches to recovering surface shape from shading, texture, 
and contours [Ikeuchi and Horn, 1981; Horn and Brooks, 1985; Render, 1980; Witkin, 1981; Brady 
and Yuille, 1984]. 

The combined output of the shape estimation processes is best collected into intermediate 
representations of the 3D shapes and configurations of visible surfaces, which we will refer to 
as visible-surface representations. Notable among proposed visible-surface representations arc the 
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depth and needle maps of Horn [1982], the intrinsic images of Barrow and Tenenbaum [1978], and 
die 2|-D sketch of Marr and Nishihara [1978]. For humans, the perception of visible surfaces is 
generally immediate, involuntary, and seems to precede (object) recognition. ITiis strongly suggests 
the existence of a visual process that autonomously computes visible-surface representations. Aside 
from the perceptual evidence, the availability of explicit visible-surface representations can also 
substantially facilitate subsequent surface analysis tasks in machine vision. 

Since early visual processing provides relative surface shape estimates with respect to the 
viewer, it is most natural to define the basic shape primitives of visible-surface representations in a 
viewer-centered coordinate system. Moreover, the primitives should be computationally compatible 
with the local depth and orientation measurements (as well as discontinuities) that are provided by 
the various shape estimation processes. These criteria are satisfied by a particularly appealing class 
of local, piecewise shape primitives known as finite elements. 

A crucial realization is that shape estimates can be provided at multiple resolutions. Indeed, 
multiresolution spatial frequency channels have been identified psychophysically in the human 
visual system (c.g., [Braddick el al , 1978]). Their existence has influenced the design of early 
visual algorithms (e.g., [Marr, 1982]). In addition, machine vision research has demonstrated that 
multiresolution processing effectively bridges fine and coarse image structure, while it simultaneously 
increases computational efficiency (e.g., [Roscnfeld, 1984]). Hence, a multiresolution organization 
of visible-surface representations is most desirable [Terzopoulos, 1982, 1983a]. 

1.2. Key Problems of Visible-Surface Reconstruction 

The main topic of concern in this paper is the development of a visible-surface reconstruction process 
responsible for generating and dynamically maintaining visible surface representations. Whether 
the intention is to model human vision or to design competent artificial vision systems, this process 
must solve four key problems [Tcrzopoulos, 1983b, 1984]: (i) the constraint integration problem, 
(ii) the interpolation problem, (iii) the discontinuity problem, and (iv) the computational efficiency 
problem. We elaborate on each of these problems next. 

(i) The Constraint Integration Problem: Each specialized visual process may be thought of as a 
quasi-independent source of information partially constraining the shapes of visible surfaces. 
The human visual system is reliable and robust because it integrates die various processes, 
enabling them to complement one another. The integration of multiple sources of information 
introduces redundancy, which is necessary not only to resolve potential ambiguities, but also 
to overcome the detrimental effects of noise and inaccuracies in the initial shape estimates. 
The constraint integration problem is fundamentally one of devising an effective means of 
integrating all available surface depth and orientation constraints (and discontinuities) within a 
cooperative visible-surface reconstruction process. 

(ii) The Inteipolation Problem: It is widely accepted that initial descriptions of images ought to 
make explicit die occurrence and local 2D structure of image features that arc correlated to 
salient events on physical surfaces (markings, boundaries, etc.). This is the essence, for instance, 
of Mart’s “primal sketch” representation of significant image irradiancc changes (edges) [Marr, 
1982]. Generally, such salient features do not occur everywhere over die visual field. The 
initial representation of images as a sparse set of features implies diat surface shape constraints 
generated by die specialized processes will also be scattered over a subset of image points. It 
is fascinating, however, that the human visual system systematically interprets visual stimuli 
such as sparse random dot stereograms as coherent 3D surfaces [Julcsz, 1971]. Indeed, diesc 
stereograms continue to elicit perceptions of dense surfaces, even when die density of dots 
carrying disparity information has been reduced until depth is unspecified over 98 percent of 
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the visible surface area (sec Fig. 1). It therefore appears that the surface reconstruction process 
is smoothly “filling in the gaps.” This phenomenon has been the subject of some psychophysical 
investigation (c.g., [Collett, 1984]). The interpolation problem of visible-surface reconstruction 
challenges us to devise a scheme, consistent with human perception, for propagating shape 
information into indeterminate regions (devoid of shape estimates) from places where it is 
available. 



Figure 1 . A sparse random dot stereogram. Binocular fusion of this stereogram reveals a planar surface as a 
central, opaque, textured square suspended nearer in depth over a similarly textured background. Vivid depth 
discontinuities separate the dense surfaces. 


(iii) The Discontinuity Problem: Visual discontinuities result from significant, spatially-localized 
changes in the physical world, particularly abrupt changes in surface structure. Both depth 
and orientation discontinuities arc perceptually relevant and provide vital boundary conditions 
for surface reconstruction. Discontinuities in depth occur at occluding contours, along which 
a surface in the scene occludes itself or another surface. Orientation discontinuities occur at 
creases or cusps of an otherwise continuous surface. In addition to the perception of coherent 
surfaces, random dot stereograms elicit vivid perceptions of surface discontinuities at abrupt 
disparity changes (see Hg. 1). The discontinuity problem amounts to (1) finding both depth 
and orientation discontinuities in surfaces, and (2) dealing with their presence during visible- 
surface reconstruction; i.e., allowing discontinuities to limit the otherwise smooth interpolation 
of shape constraints. 

(iv) The Computational Efficiency Problem: Visible-surface reconstruction at the resolution of the 
image imposes an immense computational burden on both biological and artificial vision sys¬ 
tems. Nevertheless, visible-surface representations must be computed quickly if they are to 
be of any practical value. It is generally accepted that to achieve the necessary performance, 
visual algorithms and mechanisms must emphasize parallelism [Ballard el al., 1983]; however, 
visible-surface reconstruction is compute bound to the point where the fundamental limitations 
of massively parallel mechanisms, particularly with respect to global intcrproccssor commu¬ 
nications, lead to severe inefficiencies. The computational efficiency problem is to develop a 
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visible-surface reconstruction process that not only exploits parallelism, but also overcomes co¬ 
operative communication bottlenecks to compute visible-surface representations quickly, given 
suitable architectures. For the reasons outlined in the foregoing section and following our 
previous work [Terzopoulos, 1982, 1983a], our solution to this problem hinges on the idea of 
multircsolution structuring of visual representations and associated cooperative processes. 

1.3. Prior Work 

There has been some prior work relating to the computation of visible-surface representations. 
Barrow and Tcncnbaum [1979] describe an approach to reconstructing smooth surfaces from noisy 
visual data. This approach did not apply to general classes of surfaces, however, and the proposed 
relaxation algorithms were not supported by a firm mathematical analysis. Nevertheless, Barrow 
and Tcncnbaum’s [1978] basic model of intrinsic images and much of the philosophy underlying 
their computation seems appropriate, and it has influenced our approach. 

The interpolation problem is related to classical spline approximation. A number of well-known 
surface approximation methods for scattered data are reviewed by Schumakcr [1976]. Grimson 
[1983] employed one of these methods for the continuous interpolation of visual surfaces from depth 
constraints; a minimization scheme involving a particular functional containing second derivatives 
(he referred to it as the “quadratic variation”). Brady and Horn [1983] observe that this functional 
is related to the bending energy of a thin plate (a connection noted by Duchon [1977]), and the 
thin plate model was developed further by Terzopoulos [1983a] (see also [Blake, 1984]). 

Interestingly, thin plate interpolants have appeared in other areas, including the interpolation 
of aircraft wing deflections [Harder and Dcsmarais, 1972], interpolation of meteorological fields 
[Wahba and Wendclberger, 1980], and the interpolation of digital terrain maps [Briggs, 1974; 
Bolondi el al., 1976]. In this latter paper there is some concern for the presence of discontinuities 
(faults). 

Following Ullman [1979] and others, Grimson [1983] pursued “biologically feasible,” parallel 
and iterative algorithms for surface interpolation. A serious drawback of algorithms which satisfy 
these criteria is that they often converge excruciatingly slowly for problems of reasonable size. 
The idea of multircsolution surface reconstruction exploiting multigrid relaxation methods was 
shown to overcome this problem while adhering to biological feasibility [Terzopoulos, 1982, 1983a]. 
The multiresolution methodology yields efficient algorithms not only for the surface reconstruction 
problem but for other visual problems as well [Terzopoulos, 1984]. 

In retrospect, although progress has been made, a satisfactory computational theory of visible- 
surface representations has been elusive. This is largely a consequence of the significant technical 
obstacles encountered in devising formal solutions to all four key problems of visible-surface 
reconstruction within a unified computational framework. The difficulty of the task appears to have 
evoked some skepticism as to the actual computability (hence, even the usefulness) of intrinsic 
surface representations [Witkin and' Tcncnbaum, 1983]. Based on the theoretical generality of 
our approach and the accompanying empirical results, however, we believe such skepticism to be 
premature. 


2 . 


Mathematical Analysis of Visible-Surface Reconstruction 


Let the true distance from the viewer to visible surfaces be given by the function Z(x,y), where 
x and y arc the image coordinates. Low-level visual processes generate a set of noise corrupted 
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surface shape estimates (i.c., constraints) {q} which can be expressed in the abstract notation 

c.i = £i[Z(x,y)] + (1) 

where Ci denote measurement functionals of Z(x,y) and denote associated measurement errors. 
Stated simply, the visible-surface reconstruction problem is to reconstruct, as faithfully as possible, 
the depth function Z(x,y) from the available constraints {c t }. 


2.1. The Ill-Posed Nature of the Problem and Regularization 

The problem is made nontrivial by the nature of the constraints. First, constraints are contributed not 
by one, but by multiple specialized early visual processes. Hence, slightly inconsistent measurements 
provided by different processes that happen to coincide will locally overdetermine surface shape. 
Second, constraints arc not dense, but scattered sparsely over the visual field. Therefore, while 
they may restrict surface shape locally, they do not determine it uniquely everywhere; there remain 
very many feasible surfaces that arc consistent with the constraints. Third, the measurements are 
subject to errors and noise. High spatial frequency additive noise, regardless how small its (RMS) 
amplitude, can locally perturb the surface (orientation) radically. 

In view of the above three considerations, we cannot conclude in general that the solution 
will exist, nor that it will be unique, nor that it will be stable with respect to measurement 
errors. Mathematical problems for which the existence, uniqueness, or stability of solutions cannot 
be guaranteed a priori are said to be ill-posed [Tikhonov and Arsenin, 1977], Visible-surface 
reconstruction can thus be characterized as a fundamentally ill-posed problem. 

Ill-posed reconstruction (or inverse) problems are the rule rather than the exception in early 
vision [Poggio and Torre, 1984]. Ill-posed problems cannot be solved in general, without imposing 
some additional restrictions on possible solutions. 1'his is the basis of a number of systematic 
approaches, notably the regularization methods introduced by Tikhonov and others (see [Tikhonov 
and Arsenin, 1977] and references therein). Duda and Hart [1973, Sec. 7.4] mention a basic form 
of regularization (essentially spatial smoothing) for combating the effects of noise in images. A 
more sophisticated class of regularization methods is discussed in the context of low-level vision by 
Poggio and Torre [1984]. 

Through regularization, ill-posed problems can be solved by reformulating them as variational 
principles that are effectively computable. Unlike the original problems, the variational principle 
formulations are well-posed\ i.c., it is possible to guarantee the existence, uniqueness, and stability 
of their solutions under nonrcstrictive conditions. Reformulation proceeds with the introduction of 
suitable stabilizing functionals , notably the class of stabilizers proposed by Tikhonov and Arsenin 
[1977, pp. 69-70]. These stabilizers can be interpreted as spline functionals that impose smoothness 
assumpt ions on the admissible solutions (by restricting them to Sobolev spaces of smooth functions). 1 
Pragmatically then, this type of regularization is essentially equivalent to optimal approximation by 
generalized splines [Terzopoulos, 1985a]. We pursue the generalized spline approximation point of 
view, since splines arc familiar and since they suggest helpful physical interpretations. 


2.2. A Variational Principle 

The abstract theory of optimal spline approximation is well-developed and a close connection has 
been established with variational principles involving the constrained minimization of (semi-) norms 
in (semi-) Hilbert function spaces [Laurent, 1972], Let )1 be a linear space of smooth functions and 
let $(?;) be a functional defined on V which measures the (lack of) smoothness of a function in 


1 Generic smoothness assumptions are generally the weakest (least committal) assumptions that one can make 
about feasible solutions and still obtain well-posed formulations. 
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P. Furthermore, let P be a functional on P which provides a measure of the discrepancy between 
the function and the given constraints. Consider the following variational principle: 

VP: Find u <E F such that 

£( u ) = inf <£», (2) 

where the energy functional 

C(v) - S(v) + P(v). (3) 

This variational principle will serve as a formal statement of the visible-surface reconstruction 
problem: The best reconstruction of the depth function Z(x,y) from the available constraints will 
be given by the solution u(x,y), the smoothest function in the admissible space P which is most 
compatible with the constraints. 

Before proceeding to specify the smoothness functional S (?;) and the penalty Junctional P{v), 
it should be noted that, if the solution exists, it satisfies the necessary condition for the minimum 
given by the vanishing of the first variation 6, 

6S(u) = 8S(u) + 6P(u) = 0, (4), 

which expresses the so called Euler-Lagrange equations. 

2.3. Generalized Spline Functionals 

For an appropriate smoothness functional S(w), we turn to the multidimensional splines studied 
by Duchon [1977] and Meinguet [1979], generalizations of the classical univariate splines [Ahlberg, 
et al, 1967]. The subclass of (2D) surface splines relevant to our problem can be characterized as 
members of a suitable space of admissible functions v(x,y) which minimize the functional 



The positive integer m dictates the order of the partial derivatives that occur in the functional, 
which in turn determines the order of continuity possessed by the admissible functions. The Euler- 
Lagrange equation satisfied by the minimizing function u(x,y) is an iterated version of Laplace’s 
equation: (-l) m A m w = 0, where Aw = u xx + u yy is the Laplacian of u. 

Low order surface splines have interesting physical interpretations involving equilibria of elastic 
bodies. Two special cases arc of interest. For m — 1 the functional reduces to 

Mi = / J («, + v l) dxdy, (6) 

which is proportional to the small deflection energy of a membrane (c.g., rubber sheet), while for 

m = 2, 

Hi - j J (v 2 xx + 2 v\ y +'«*„) dxdy, (7) 

is proportional to the small deflection bending energy of a thin plate (with zero Poisson ratio) 
[Courant and Hilbert, 1953]. Duchon [1977] refers to the minimizers of \v\? 2 as thin plate splines. 

Since thin plate splines arc the natural 2D analogs of cubic splines, \v\\ finds frequent usage in 
surface interpolation problems [Schumaker, 1976]. In particular, it has been employed for visual 
surface interpolation [Grimson, 1983; Terzopoulos, 1983a]. 

The physical interpretations make it clear that membrane splines offer a lower order of 
continuity than thin plate splines. Since the physical forces in the membrane arc due primarily 
to its surface tension, it generates minimal area surfaces. Although minimal area surfaces are 
continuous, they need not have continuous first partial derivatives; i.c., they are C° surfaces. For 
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instance, a sharp corner would result readily if an idealized physical membrane were subjected to 
the deflecting force of a knife edge. In contrast, the restoring forces in a physical thin plate arc 
due primarily to its flexural rigidity. A thin plate would not crease when deflected by a knife edge. 
Thin plate splines therefore maintain continuity as well as continuous first partial derivatives; i.e., 
they generate C l surfaces. 


2.4. Controlled Continuity and the Thin Plate Surface Under Tension 


Generic smoothness assumptions arc justified in pursuing a regularization approach to the visible- 
surface reconstruction problem, inasmuch as the coherence of matter tends to give rise to smoothly 
varying surfaces relative to the viewing distance, over some range of scales; however, smoothness 
assumptions clearly do not hold arbitrarily across surface discontinuities, some of which persist 
across all scales. This introduces significant complications for classical spline approximation or 
regularization methods; the continuity of spline functionals (or stabilizers) must be controlled at 
discontinuities in order to preserve them. 


A stabilizer providing the necessary local continuity control can be realized as a weighted 
combination of generalized spline functionals of more than one order m [Terzopoulos, 1985a]. We 
propose the following smoothness functional; 


S pT {v) = j J J^p(x,y)[r(x,y)(vl x + 2v 2 xy +u“ y ) + [1 - t(x, y)](v 2 x + vj|)} dx dy, (8) 

where 0 denotes the image domain, and p{x,y) and r(x,y) are real-valued weighting functions 
whose range is [0,1). This controlled-continuity stabilizer is a weighted convex combination of the 
thin plate spline functional |u |2 and membrane spline functional |v|J integrands. The associated 
Euler-Lagrange equation is 


d 


dx 2 


, \ & / n \ & i \ & t \ 

( P'Uxx) + oftf— {2P' u xy) + £^2 \^ U yy) ~ fff ( T l u x) 


d_ 

dy 


(rjUy) = 0 , 


( 9 ) 


where p(x,y) = p(x,y)r(x,y) and rj{x,y) = p(x,y)[l-T(x,y)}, with natural (i.e., free) boundary 
conditions. The functional 5 / , 7 -(f) can be thought of as a thin plate surface under tension , where 
p[x,y) is a spatially varying “rigidity” and [1 - r(x,y)\ is the spatially varying “surface tension.” 
It generalizes the unidimensional splines under tension of Schwcikert (see [Ahlberg, et ah , 1967]). 


The local continuity properties of the thin plate surface under tension functional can be 
controlled at any point (x,y) e 0 by specifying the values of the continuity control functions 
p(x,y) and r(x,y) at that point. As t approaches 1 the functional tends to a thin plate spline (a 
C l surface) whereas towards the other extreme, 0, the functional tends to a membrane spline (a 
C° surface) with intermediate values characterizing a hybrid C 1 surface that blends the properties 
of both constituent splines, p determines the overall potency of the smoothness functional. 


Reconstructed surfaces must be able to faithfully preserve known depth and orientation 
discontinuities, while not introducing spurious discontinuities at other locations. This can be 
accomplished if (i) away from known depth'and orientation discontinuities, the reconstructed 
surface possesses (at least) the C 1 smoothness of a thin plate, maintaining both continuity and 
continuous derivatives, (ii) at known orientation discontinuities, it exhibits just the C- 1 smoothness of 
a membrane, maintaining continuity only, and (iii) at known depth discontinuities, the smoothness 
functional is deactivated so that the reconstructed surface is free to “fracture” locally. Hence, 
$ pT (v) will be manipulated as follows: At all non-discontinuity points (x,?/), p{x,y) and r(x,y) 
should be nonzero. At orientation discontinuity points, r(x,y) is set to zero. At depth discontinuity 
points, p{x,y) is set to zero. Mechanisms for automatically detecting discontinuities by computing 
continuity control functions optimally according to local criteria are considered in a subsequent 
section. 




Terzopoulos 


Computing Visible-Surface Representations 


8 


2.5. Penalty Functionals 


Assuming independently distributed measurement errors c t with zero means and variances af, the 
optimal measure of incompatibility is a weighted Euclidean norm of the discrepancy between the 
admissible function and the data c t -: 

( 10 ) 

i 

where the a,- are nonnegativc real-valued weights (ideally a, is inversely proportional to of; 
i.e., = 1/A of) [Kimcldorf and Wahba, 1970], This penalty functional can also be employed 
(suboptimally) when the above assumptions do not hold strictly. 


Appropriate measurement functionals Li for surface reconstruction may be synthesized from 
generalized /c^-ordcr derivatives: 



d k v 


dx3dy k ~i 


1 



(ii) 


k = 0 yields simple evaluation functionals £i{v(x,y)] = v(xi,yi), which will be employed to 
model the local depth constraints 

c t = v(xi,yi) + ci = d {x . iy .y (12) 

The components of the local surface normal n [x{,yi) = [v x (xi,yi),v y (xi,yi), -1], which deter¬ 
mine local surface orientation, can be handled by the first order (k = 1) derivative functionals 
£ t -[o(x,y)] = v s (x,-, yi) and Ci[v(x,y)] — v y (xi ,?/ t ) and yield analogous expressions for the local 
orientation constraints: 

Cj =»,(!<,») + «.'=?(*„»<) 

^ v y(x,.y,) 

Other potentially relevant functionals such as directional derivatives can be accommodated straight¬ 
forwardly with the above notation. 


It is convenient to separate the various constraints into three sets; the set i £ D of image 
points at which depth constraints occur, and the sets i £ P and i £ Q at which orientation 

constraints P( Xx , y ,) and q( Xi , yi ) occur respectively. The penalty functional can then be expressed as 
a sum of three components 




ieo 


i€P 


i£Q 


where the o t - parameters are now distinguished as n Pi> and a q{ . 


(14) 


2.6. A Physical Model for Visible-Surface Reconstruction 

The variational principle formulation of the surface reconstruction problem has an appealing 
physical interpretation which is illustrated in Fig. 2. The thin plate surface under tension may 
be visualized as an elastic surface, planar in its natural state, whose elastic bending energy S pT (v) 
stabilizes surface shape so that it varies smoothly in between constraints (but not at discontinuities). 
Constraints deflect the surface according the penalty functional P(v), which can be interpreted as 
the total stretching energy of a set of ideal springs attached to the constraints. The left part of 
the figure shows the clastic surface whose deflection u(x,y) at equilibrium is determined by an 
infrastructure of scattered depth constraints. The local depth estimate is encoded as the vertical 
height of the constraint and the tightness of each constraint is controlled by associated spring 
stiffness a ( ii- The right part of the figure illustrates an orientation constraint coercing the local 
surface normal. The spring stiffness is determined by the constraint parameters a Pi and a qi . 
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Figure 2. The physical model. Thin plate surface under tension and depth constraints (left). Local influence 
of an orientation constraint (right). 


2.7. Existence, Uniqueness, and Stability of the Solution 

Existence, uniqueness, and stability of the solution u(x, y) to the variational principle VP are 
guaranteed when £ pr (v) = Spr(v) + P(v) is a norm in the admissible space M. Unfortunately, 

A 

generalized spline functionals |u| m are a priori only semi-norms (of a particular class of Sobolev 
spaces). The null spaces M of functions that map to zero under the semi-norms are simply the 
(M = ("V -1 ) dimensional) spaces of all polynomials over 9? 2 of degree less than or equal to m --1 
[Meinguct, 1979], The penalty functional P(v) can force £ pr (v) to be a norm, however, if it at 
least constrains M to a unique polynomial. A possible set of conditions for this to occur is that the 
Li include evaluation functionals at an M- -unisolvent set of points (i.e., a set of M points which 
define a unique polynomial in the null space of the smoothness functional). In particular, since the 
maximum order of generalized splines in the stabilizer S pr (v) is m - 2, its null space is the space 
of linear polynomials. Thus, the following proposition can be proven [Terzopoulos, 1984]: 

Proposition. The solution u(x, y) will exist, be unique, and stable given any one of the following 
minimal conditions 

(i) three noncolinear depth constraints, 

(ii) two depth constraints as well as a single p or q constraint, 

(Hi) a single depth constraint as well as a single.p and a single q constraint, 

(iv) a single p and a single q constraint with the “center of gravity ” of the surface fixed. 

These minimal conditions will hold in practice, due to the large number of constraints typically 
available from early shape estimation processes (the fixed center of gravity condition can be imposed 
when necessary). Consequently, the visible-surface reconstruction problem may be considered well- 
posed, hence effectively computable in general. 

Satisfying the conditions for a well-posed problem essentially guarantees that a unique state of 
stable equilibrium will exist for the platc/spring system (the minimal energy shite £ pr (u)). In this 
context, the controlled continuity assumption about surfaces, as embodied by the thin plate surface 
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under tension model, is physically nonrcstrictivc but nonetheless powerful enough to guarantee the 
existence of unique solutions to die variational principle. 


3. Discretization 


It is extremely difficult, if not impossible, to obtain an analytic solution to the variational principle 
due to the irregular occurrence and geometry of constraints and discontinuities. For our purposes, 
the only viable approach is to convert the continuous surface reconstruction problem to an equivalent 
discrete problem whose solution can be computed numerically. To this end, finite elements make 
ideal local surface shape primitives for use in visible-surface representations [Terzopoulos, 1982, 
1983a]. The finite element method [Strang and Fix, 1973] is a general, powerful, and mathematically 
rigorous approximation technique which guides the selection of appropriate elements and governs 
their interactions according to the nature of the variational principle. 


The finite element method offers substantial flexibility in discretizing domains with irregular 
shaped boundaries. Although the use of irregularly shaped elements to discretize such domains 
may not present a feasibility problem with regard to distributed biological mechanisms, it makes 
nontrivial the mapping of elemental computations onto regularly interconnected processing networks 
typically provided by VLSI technology. In this paper we restrict ourselves to regular finite elements 
in order to facilitate such mappings. Since the goal is to obtain a particularly fine discretization, at 
the resolution of the image, the restriction to fine regular elements will not jeopardize our ability 
to accommodate the irregular occurrence of constraints or discontinuities. 


3.1. The Discrete Equations 


The domain is tessellated into square element subdomains with sides of length h. Nodes are 
located at clement corners and shared by adjacent elements. This results in a planar and uniform 
square grid of nodes that is ideally suited to VLSI implementation. The nodes arc naturally 
indexed by (*, j) for i = 1,..., TV* and j = 1,..., N y , where N x and N y are the number of 
nodes along the x and y axis respectively of the (rectangular) domain 0. The total number of nodes 
is N - N x x N y . The reconstructed surface is represented by an assembly of (nonconforming) 
finite elements, each of which is a six-point (full) quadratic intcrpolant defined locally within its 
particular subdomain. The unknown displacement (surface depth), at node (i, j) is denoted by the 
variable v£ ■ = v(ih,jh). Taken together, the displacement variables arc denoted by the vector 

\ h e Once this vector is determined by solving a discrete version of the variational principle, 
the local intcrpolants arc known exactly and, consequently, they explicitly represent depth and 
orientation everywhere over the surface. 


The proposed square, quadratic clement leads to the following 0(/r) formulas for the required 
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partial derivatives at an arbitrary node (i,j) [Terzopoulos, 1983a]: 


^ h 2 


v 


yy ^ h 2 


^ (v?, 1.3 “ 2v ?,3 + v r 1 , 3 ) ; 
4(vf,3 (1 -2v?, 3 + v^ 1 ); 


zy j t 2 ( v t-t 1,3 + 1 V i,3 + 1 V * + l,3 V +j) 1 

h ( V, * +1 J ” dl) ’ 


(15) 


17 ^ =£■ 

* h 


vh 
» h 


Note tliat the formulas arc finite difference expressions. Their appearance is due to the uniformity 
and low order of the clement and their relative simplicity will facilitate the calculations substantially. 

Substituting the above expressions with the constant approximations p(x,y) => p^ ■ and 

r(x,y) => r£. into (8), and noting that the area of each clement is h 2 , we obtain the discrete 


functional 




JCl 

h 2 




+ 


1 - T: • 


- h ] ( 

h3 \ 


( V ?+l,3 “ 2V ?,3 + Clf 2 

+ 2 (vJ+i.j+i - v?j- +1 - v? +1>J . + *?„■) 
+ (v?, 3 +i - 2v? >y + V^ . ,) 

\i) 2 + ( V ?,3+l - V ?,3') 


(16) 


V ?+l,y “ V i,3 


Although by no means a necessity, it is both natural (in view of common image discretization) 
and convenient to assume that the constraints coincide with nodes (i,j) of the grid. Hence, to 
obtain a discrete expression for P(v ), we collect the nodes at which the various constraints occur 
into three sets; die set (i. j) £ D at which depth constraints dJ* • occur, and the sets (i,j) £ P 

and (i,j) £ Q at which orientation constraints p£ • and q|* • occur. Using symmetric difference 
approximations for the partial derivatives in (13), the discrete penalty functional may be written in 
terms of the nodal variables as 



The energy-minimizing vector of nodal displacements u h satisfies the equilibrium condition 

V#(u fc ) = VS f ; T (u h ) + WP h (u h ) = 0, (18) 

where V is the gradient operator. Since the discrete functional Sp T { u /l ) is a quadratic form in the 

uf -, the above equation defines a linear system of simultaneous equations that are satisfied by u h . 
The discrete problem amounts to solving these nodal equations. 
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3.2. Computational Molecules 


To progress towards explicit expressions for the nodal equations, we first determine the partial 
derivatives of Sp T (u h ) and P h (u h ) with respect to an arbitrary nodal variable ufy. Letting 

(19) 


h _ Ji _ft i -i „ft _ _ft (t _ft 


we obtain 


Hi = 


and 


rt h . 

Ii,3 


P T ' 

du^ ■ 


u ?,y “ 2u ?-i,y + * 


+ 

+ 

+ 


t+lj 


+ u t,y) ^i+i,y 


<-*,)) ri-lj 

{~ 2 a i+i ,y + 4u ?,y - 2u ?-i,y) ft\y 

( u ?+2,y - 

( 2u L - 2u .--i,j - 2 u L-i + 2u ?-i,i-i) 

+ (- 2u ?+i,j + 2u ?,y + 2 «?+i,y-i - 2 <y-i) fty-1 
+ (- 2 u L+i + 2u ?-i,y+i + 2u ?,y - 2u ?-i,y) /‘.' lj 
+ ( 2u ?+i,y+i - 2u ?,y+i - 2u ?+i,y + 2«?,y) 

+ ("t ~ 2 <y-i + <y-j) <-i 
+ (~ 2u .\y+i + 4u ?,y - 2u ?,y-i) ft-.y 
+ ( u ?,y+2 - 2u .\y+i + u !j) | 

"*"( ( u i,y ~ u i-i,y)l/i-iy + ( u i,y ~ u i+i,y) ^i,y 
+ ( U L “ u ?,y-i) ’li-j-i + ( u ?,y - <y+i) >&}> 


the discrete version of (9). Next, for (i,j) e D n P ft Q, 


, / a Pi~U (..h ft \ _ ft 

+ [ 4 h 2 \ i ~ 2 ^) 2 h 


4. ( + ( u h ft \, a PiHd h 

^ 4 h 2 \ '- J + 2 h Pi+1 ’* 

+ (~4^~ K - U ‘^ 2 ) - -2T- q ^- 1 

/ „ ft ft 

, / ghl±l /..fc fa 

+ \ 4/i 2 \ l ’ 5 ' + 2 h 


( 20 ) 


( 21 ) 


The above expressions specify the nodal equations implicitly. Each constituent term in (round) 
parentheses can be represented graphically as a basic computational molecule. Computational 
molecules will be interpreted both as spatial representations of the nonzero coefficients in die nodal 
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equations, and as local computations involving multiplications and additions of specific proximal 
nodal variables. The former interpretation will facilitate the construction of individual nodal 
equations given some local structure of constraints and discontinuities, while the latter will lead 
directly to local iterative algorithms for solving the resulting simultaneous linear system. 


3.2.1. Basic Molecules 


A 



Figure 3. Plate molecules (a) and membrane molecules (b). 
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Kq. (20) is a convex combination of two components; the first stemming from the thin plate 
energy functional, the second, from the membrane energy functional. Kach constituent term yields 
a basic computational molecule (sec Fig. 3). a set of linked atoms indicated by circles. The central 
node (i,j) is indicated by a double circle in each molecule. Fig. 3(a) illustrates the ten plate 
molecules obtained from the terms of the first component, while Fig. 3(b) shows the four membrane 
molecules obtained from the terms of the second component. Each atom contains the coefficient of 
the associated nodal variable (aside from die /if. and r?f. factors). 


A 



B 





Figure 4. Depth constraint molecule (a) and orientation constraint molecules (b). 

Similarly, the depth constraint term in (21) can be represented by the depth constraint molecule 
shown in Fig. 4(a). Associated with it is the factor oyfydf y, which is indicated underneath the 
molecule. The remaining orientation constraint terms of (21) arc represented by the orientation 
constraint molecules and associated factors shown in Fig. 4(b). 
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3.2.2. Molecular Summation within Smooth Regions 

The formation of nodal equations within continuous regions can be visualized as a process of 
molecular summation. During molecular summation, the basic molecules combine at the central 
node, coincident atoms summing together. 

When (t,j) is an interior node, away from constraints and discontinuities, p, h - = r^. = 1, and 
only the plate component of (20) contributes to the expression for the partial derivative. Hence, 
the equilibrium condition (18) reduces to the nodal equation 

0 = j? u L ~ £* ( u ?-u + u !Vu + 1 + 

+ p ( u ?-i,j-i+ u ?+i,y-i + u i-i,j+i + u !+i,y+i) ( 22 ) 

+ p ( U ?-2,J + U i+2,j + u i\y-2 + ^+ 2 ) • 

This equation can be represented by the composite nodal molecule illustrated in Fig. 5(a), which 
results from the summation of the plate molecules in Fig. 3. 




A 


Figure 5. Interior node molecules, (a) Away from discontinuities, (b) At interior orientation discontinuities. 

Note that the computational molecule for the center of the region is a factor of h 2 (due to die 
elemental area) times an order 0(/t 2 ) finite difference approximation for die biharmonic operator 
[Abramowitz and Stegun, 1965, p. 885], the Eulcr-Lagrange equation associated with the thin plate 
spline. This is an expected consequence of die particular clement employed which yielded finite 
difference approximations for the second partial derivatives of v h . 

If node (i,j) is a depth constraint, the first term in (21) takes part in the nodal equation. The 
effect can be represented as a summation of the depdi constraint molecule and associated constraint 
factor with the nodal molecule for (i,j) shown in Fig. 5(a). 

Similarly, if (i - 1 ,j) or (t +1 ,j) arc p constraints, or (i,j - 1) or (*, j T1) arc q constraints, 
the other terms in (21) participate in the nodal equation. Again dais can be represented as the 
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summation of computational molecules. Specifically, the upper left molecule in Fig. 4(b) sums with 
the nodal molecule if (i - 1 J) e P. the upper right only if (i + 1 ,j) e P , the lower left only if 
[i,j - 1 ) eQ, and the lower right only if (i,j + 1 ) e Q. 

3.2.3. Molecular Inhibition at Discontinuities 


If (t, j) is a discontinuity, cither pf>- or 17 J . or both may be zero, thus nullifying the summation 
of specific molecules. This crucial influence of the discrete continuity control functions near known 
discontinuities will be referred to as molecular inhibition. It was convenient for expressing (20) and 
(21) to discretize u(x,y), p(x,y), and r(x,y) over the same set of nodes. Although orientation 
discontinuities can be situated at these nodes (since uf • is defined at an orientation discontinuity), 

it is better to position depth discontinuities on the links half way between nodes (since u* ■ is 
undefined at a depth discontinuity). Discretizing t (x, y) on links does not present a problem in 
practice. As a general rule, a discontinuity may inhibit a molecule only if it coincides with a 
constituent atom or link. 


First consider orientation discontinuities. At an orientation discontinuity, r*. = 0 and only 
the second component of (20) contributes to the nodal equation. In effect, the plate molecules are 
inhibited and replaced by the membrane molecules of Fig. 3. At an interior orientation discontinuity 
away from depth discontinuities, all four membrane molecules superpose to yield the nodal 
molecule shown in Fig. 5(b), which represents the nodal equation 


4u,- 


M 


u- 




n h 

u t+i,J 


u, 


»\y-i 


- u 


«,y+i 


0 . 


(23) 


The equation will be recognized as - h 2 times a standard finite difference equation for the Laplacian 
[Abramowitz and Stegun, 1965, p. 885]. It too appears as a consequence of the Euler-Lagrange 
equation associated with the membrane spline. 


Since an orientation constrain cannot meaningfully coincide with an orientation discontinuity, 
orientation discontinuity nodes inhibit orientation constraint molecules. On the other hand, depth 
constraint molecules are not inhibited by orientation discontinuities since it is perfectly reasonable 
to locally constrain a membrane spline in depth. 

Because smoothness constraints arc unsuitable at a depth discontinuity node (i,j) (i.e., pij = 
0 ), a nodal equation cannot involve nodal variables separated by or coinciding with a depth 
discontinuity. Consequently, depth discontinuities inhibit all of the basic computational molecules. 
Fig. 6 (a) illustrates examples (disregarding constraints) of nodal molecules for boundary nodes 
(marked as double circles) which are near depth discontinuity nodes (marked by X’s). Examples of 
nodal molecules at boundary orientation discontinuities (double circles) next to depth discontinuities 
(X’s) are shown in Fig. 6 (b). 


4. Detection and Localization of Surface Discontinuities 


An important feature of our framework for computing visible-surface representations is the uni¬ 
form treatment of constraints and discontinuities, essentially as localized and independent surface 
shape primitives. This facilitates the parallel integration of discontinuity information, along with 
shape constraints, over the various early shape estimation processes. It is convenient to think of 
discontinuity information as being collected into a discontinuity map which is in registration with 
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Figure 6. Molecular inhibition at discontinuities, (a) Nodal molecules at boundary nodes (double circles) near 
depth discontinuity nodes (X's). (b) Nodal molecules at boundary orientation discontinuities (double circles) 
next to depth discontinuities (X’s). 
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the reconstructed surface. Technically, the map comprises the nodal variables {/??•} and {r.M 
representing the discrete continuity control functions. 

Any early visual process can participate in initializing die discontinuity map according to its 
own local hypotheses about the occurrence of discontinuities. In general, this prior discontinuity 
information will be partly incomplete and inconsistent, since it derives from narrowly specialized 2D 
image analysis. In evolving a globally consistent surface, the visible-surface reconstruction process 
performs a crucial task: it brings the prior discontinuity information into consonance with the 3D 
shape constraints collected from all the early processes. 1'his raises the problem of detecting and 
localizing both depth and orientation discontinuities during reconstruction. The current section 
investigates this problem for the (impoverished) case in which no prior discontinuity information 
is available. We first propose a straightforward scheme which exploits the regularizing properties 
of the surface model, then a more sophisticated approach that extends the variational principle to 
optimally estimate discontinuities according to generic expectations about their local structure. 

4.1. Regularization Based Discontinuity Detection 


From one perspective, surface discontinuity detection shares much in common with traditional 
approaches to image intensity edge detection. In particular, it is possible to detect discontinuities 
by applying thresholded local differencing operations to the reconstructed surface which, like the 
image, is a regularly sampled function. Because they are easily corrupted by image noise, however, 
local edge operators such as Laplacians perform poorly [Rosenfeld and Kak, 1982] without a 
smoothing prefilter, say a Gaussian [Marr and Hildreth, 1980]. Interestingly, the thin plate surface 
under tension performs the necessary smoothing on the sparse and noisy shape constraints (standard 
low-pass filters such as Gaussians are inapplicable to sparse data). This regularizing effect permits 
the reliable computation of numerical derivatives for detecting discontinuities [Bakhvalov, 1977, 
Sec. 5.4; Poggio and Torre, 1984; Terzopoulos, 1985a]. In addition to exploiting the regularizing 
effect of the thin plate surface under tension, die discontinuity detection scheme described next 
is easily accommodated within the distributed computational structure of our framework, and it 
permits relevant criteria such as psychophysically measured limits on stereofusion to impact on 
discontinuity detection. 

Consider the random dot stereogram in P'ig. 7 which depicts a set of planar surfaces stacked in 
depth. Fig. 8 shows a single continuous surface generated by the surface reconstruction algorithm 
from sparse stereoscopic disparities provided the Marr-Poggio-Grimson (MPG) stereo algorithm 
[Grimson, 1985]. Fig. 9 dramatizes a portion of the reconstructed surface in cross section as it passes 
across a depth discontinuity. The C 1 surface overshoots constraints near the discontinuity because 
its smoothness conflicts with the sudden change in depth. The surface is clearly inappropriate as a 
final solution near depth discontinuities, but the local incompatibility can signal the occurrence of 
these discontinuities. 


Opposing bending moments are imparted to the surface by the constraints on either side of 
the discontinuity. The surface inflection (sec Fig. 9), where the bending moment undergoes a sign 
change, localizes the depth discontinuity. For a thin plate spline u(x,y), the bending moment per 
unit length parallel to the x-z plane is proportional to while its counterpart parallel to the 
y-z plane is proportional to —u yy [Szilard, 1974]. The sum of orthogonal bending moments gives 
the total moment M = -{u xx + u yy ) = -Aw, the negative Laplacian of the deflection function. It 
can be computed readily at a node (i,j) of the discrete surface using the standard approximation: 
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(24) 



Terzopoulos 


Computing Visible-Surface Representations 


19 



Figure 7. Synthesized random dot stereogram. When fused, the stereogram depicts four planar surfaces 
stacked one atop the other in depth. 


The zero crossings of M for the reconstructed surface in Fig. 8 are shown on the left in 
Fig. 10 as black contours. Most of these correspond to weak inflections due to slight ripples in the 
reconstructed surface. A measure of significance is therefore needed to detect true discontinuities 
while weeding out spurious, weak inflections. The magnitude of the local depth gradient (surface 
tilt) is a suitable significance measure for depth discontinuities. Hence, an inflection point will 

be considered significant if G = |Vu| = ^ + u 2 exceeds a limit td (it is more efficient to 

use the square of this expression u 2 x + u 2 , or even \u x \ + |ti y |). Employing the usual discrete 
approximations, we obtain 



The right half of Fig. 10 shows the significant inflection points where Gij > td with td = 1. 
Adding these significant points to the discontinuity map (by setting the associated p 1 } - to zero) 
fractures the continuous surface to yield as a solution the reconstructed stack of surfaces shown in 
Fig. 11. 

The limit t d must be large enough so that weak inflection points are rejected as possible 
discontinuities, while not so large as to miss many true depth discontinuities. A possible criterion 
for choosing td in applications to stercopsis of opaque surfaces is suggested by Panum’s limiting 
case: i.e., when a surface is tilted so much from the viewer that it begins to occlude itself from one 
eye, causing stercopsis to fail. Human stereofusion limits have been measured psychophysically. 
Using pairs of points at different orientations, Burt and Julcsz [1980J measured a roughly isotropic 
disparity gradient limit of approximately 1 between fusion and diplopia. Interestingly, this is only 
half the Panum limit. 

It is not inconsistent with these findings to use the disparity gradient limit td to detect 
significant depth discontinuities in conjunction with the isotropic bending moments M f ;j to localize 
these discontinuities. The required local support computations can be performed in parallel at each 
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Figure 8. Single surface reconstructed from the stereogram of Fig. 7. 


grid point over the surface. Analogously, significant orientation discontinuities may be detected 
when the magnitude of the bending moment \M{j\ of the surface exceeds a limit t Q (points of 
high curvature), and they may be localized at relative extrema of the bending moment (positions 
of locally highest curvature). The sign of a bending moment extremum indicates die sense of the 
orientation discontinuity; negative signals a concave crease, and positive, a convex crease. Curvature 
peaks were also employed in a scheme for detecting surface orientation discontinuities proposed by 
Langridge [1984]. 

4.2. Discontinuity Detection by Variational Continuity Control 

On the one hand, experimentation on natural data with the regularized approach to discontinuity 
detection demonstrates the feasibility of discovering many of the more significant discontinuities 
during surface reconstruction (results are presented later). On the other hand, certain inherent 
inadequacies of this simple scheme can often lead to poor surface reconstructions. The shoitcoming 
are due to a basic conflict caused by smoothing. While regularization eliminates noise, making 
reasonable estimation of surface derivatives possible in continuous regions, it tends to obscure 
discontinuities [Tcrzopoulos, 1985a], It can result in poor detectability and localization of the more 
subtle discontinuities, a common problem with smoothing edge operators in general [Leclcrc and 
Zucker, 1984]. 






Terzopoulos 


Computing Visible-Surface Representations 


21 


l^ significant inflection 

V insignificant inflections 



Figure 9. Cross-section of a reconstructed surface across a depth discontinuity. The significant and insignificant 
inflections of the surface are indicated. 


The problem can be resolved by exploiting more fully the controlled continuity model to 
preserve surface discontinuities and, moreover, to incorporate a priori expectations about dis¬ 
continuity structure into the variational principle for surface reconstruction. So augmented, the 
variational principle establishes a beneficial cooperation between the interpolation process, which 
smoothly propagates shape information across regions, and the complementary discontinuity pro¬ 
cess, which delimits these regions. Thus it optimally reconstructs the piecewise continuous surfaces 
and discontinuities simultaneously to achieve the best possible surface shape. 

As was mentioned in the previous section, the smoothness of the thin plate surface under 
tension is incompatible with any sudden transitions imposed by the scattered shape constraints. 
This implies that its potential energy of deformation is generally greater at what ought to be 
interpreted as surface discontinuities. Any local reduction in the continuity of the surface reduces 
the incompatibility and locally reduces potential energy. This can be seen from (8); S pT (v) 
considered as a function of (v, p , r), decreases as either p(x, y ) or r(re, y ) arc made zero over more 
of Cl. This suggests that discontinuities can be discovered in die course of solving the variational 
principle, by allowing the surface to crease and fracture as needed to reduce the total energy 
below the minimum obtainable with a single smooth surface. The insertion of discontinuities must, 
however, incur some energy increase, otherwise p(x , y) — 0 everywhere would trivially minimize 
the energy. 
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Figure 10. Bending moment zero crossings (left) and detected depth discontinuities (right) for the reconstructed 
surface in Fig. 8. 


The variational continuity control approach to detecting discontinuities involves augmenting 
the original energy functional £ pT [v) = S pT (v) + P[v) with the discontinuity functional P(p,r) 
(explained shortly) to obtain the new variational principle 

Find u, p, and r such that 

6(u,p,f) = inf e(v,p,T) t (26) 

V,P,T 

where (he energy Junctional 

£{v,p,r) = S(v,p,r) + P(v) -I- P(p,r). (27) 


The solutions u(x, y), p(x,y), and r(x,y), satisfy the three coupled Eulcr-Lagrange equations, 
which express the vanishing of the first variation with respect to each independent function 

3 (^y); 


d d d d 

6 u £{u,pj) = 0 fin**) + {Zp*xy) + ^2 (P u vv) ~ ffx dy 


8 P S (u, p, r) = 0 =r(vl x + 2v\ y + vj y ) + [1 - r](u* + v 2 ) + 8 p P(p , r); 
8 r £(u,p,r) = 0 =p^(v 2 xx + 2v xy + v 2 y ) - (v\ + v 2 y ) + 8 T P(p,f). 

Note that the first equation is identical to (9). 


(28) 


The functional P(p,r) maps the depth and orientation discontinuity configurations p(x,y) 
and r[x,y) into positive energies (this is analogous to the role of S pT (v) with respect to u ). In its 
simplest form, the functional can increase monotonically with the total number of discontinuities; 
e.g., D(p, t) = ff n Ai[l “ p( x , y)] '+ A>[1 - r(x, y)] dxdy, where fi d and (3 0 are positive energy 
scaling parameters for the depth and orientation discontinuity contributions, respectively. 


More interestingly, significant advantages accrue in the detection of weak or subtle disconti¬ 
nuities if the functional can be designed so as to bias the solution according to generic constraints 
about the local structure of discontinuities. Useful constraints can, for example, be based on 
Gestalt principles of good continuation — discontinuities tend to be arranged along contours, these 
contours tend to be continuous, etc. This may be accomplished readily by assigning potential 









Terzopoulos 


Computing Visible-Surface Representations 


23 



Figure 11. Reconstructed surfaces and discontinuities. 


energies to various local discontinuity configurations on the (i,j) grid of nodes for die discrete 
problem. Encodings of local edge configurations that favor good continuation have been employed, 
for instance, in relaxation labeling curve enhancement processes [Zucker el al ., 1977] and in Markov 
random field image restoration models [Geman and Gernan, 1985]. 

In our current implementation, the discrete discontinuity functional is a weighted nodal sum 
of potential energy quanta Df- and • over depdi and orientation discontinuity configurations 
respectively: 

P V, r h ) = (P h ) + /*Ofc(r*)]. (29) 

We employ a heuristic encoding which favors die formation of continuous and smoothly curving 
contours by locally assigning higher energies to isolated discontinuities, terminations, sharp bends, 
junctions, and regions. Fig. 12 illustrates some of the configurations, and the (numeric) energies 
associated with them and their rotationaily symmetric counterparts. Just as for computational 
molecules, the circles denote nodes while die X’s denote discontinuities (positions where p h 
or r h are 0). The quanta in Fig. 12(a) constitute D^ - for depth discontinuities, which occur on 
links between nodes as explained previously. They are equivalent to die configurations found in 
[Geman and Geman, 1985]. Fig. 12(b) depicts some of die orientation discontinuity configurations 
encoded by Oj* ■. Orientation discontinuities coincide with nodes. 







The discrete variational principle simultaneously governs the values of the displacement nodal 
variables of the surface as well as die nodal variables in die discontinuity map. Although die energy 
functional S p , T (v) has a unique minimum (given die conditions of Sec. 2.7) for fixed p and r, this 
is no longer die case for <f (<;,/>, r) which allows variation of die continuity control functions in the 
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minimization. The nonconvcxity of the energy landscape makes this a much more difficult problem 
to solve numerically. In the discontinuity detection experiments to be presented in Sec. 6.5, we 
propose a strategy for efficiently obtaining good, though not necessarily optimal solutions. 


5. Overview of the Multiresolution Surface Reconstruction Algorithm 


Application of the finite element yields the discrete problem of solving a linear system of simul¬ 
taneous equations. This system has computationally desirable properties; i.e., its matrix is sparse, 
banded, symmetric, as well as positive definite (for fixed p(x,y) and tau(x,y)) when the available 
constraints satisfy the conditions for a well-posed problem. The sparseness of the matrix, a direct 
consequence of the local support of the finite element, is evident from the nodal equation of an 
interior node; rows associated with interior nodes have only 13 nonzero entries, while nodes at and 
near discontinuities have even fewer. The N x N matrix however, tends to be extremely large in 
practice, since the number of pixels N in a typical image can range from 10 4 to 10 6 or greater. 
This combination of properties suggests the application of iterative techniques such as (parallel) 
Jacobi or (sequential) Gauss-Seidel relaxation methods [Hagcman and Young, 1981]. Relaxation 
methods lead to distributed algorithms, and the parallel variants may be implemented concurrently 
on networks of many simple, locally-interconnected processors. 


5.1. Nodal Relaxation Computations 


A local-support nodal relaxation compulation can be obtained at node (t, j) by expressing u* • 
in terms of the remaining variables in the nodal equation determined by the local structure of 
constraints and discontinuities. The nodal relaxation computation may be constructed automatically 
by applying our simple rules governing the summation of basic computational molecules: 

(i) Plate, depth constraint, and orientation constraint molecules sum at interior (non-discontinuity) 
nodes. 


(ii) Membrane and depth constraint molecules sum at orientation discontinuity nodes. 

(iii) Orientation discontinuities inhibit plate and orientation constraint molecules. 

(iv) Depth discontinuities inhibit all basic molecules. 

For instance, at a depth constraint node away from discontinuities, the Gauss-Seidel relaxation 
computation becomes 
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(30) 


where we have suppressed the discretization superscript h and instead introduced the bracketed 
iteration indices. At an unconstrained depth discontinuity, we obtain 

(nfl) _ 1 ( (n + l) , (n) , „( ri + 1 ) , ,.(«) \ (oi) 

u » J - 4 \ u i 1J + 1,3 1 U t,i-1 + u i,jn) • ( 31 ) 

Note tliat the nodal relaxation computations do not change from one iteration to the next, so long 
as the influencing constraints or discontinuities remain unperturbed. 
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5.2. Multiresolution Relaxation 

A serious problem with iterative techniques, in general, is their slow convergence rates for large 
problems. This inherent inefficiency is due to the fact that information must propagate incrementally 
across large representations from nodes to their near neighbors in accordance with the nodal 
relaxation formulas. 2 We have developed highly efficient iterative algorithms that overcome this 
problem for surface reconstruction [Terzopoulos, 1982, 1983a] as well as for certain other visual 
problems [Terzopoulos, 1984]. These algorithms achieve efficiency by exploiting multiresolution 
relaxation methods [Fedorenko, 1961; Brandt, 1977; Hackbusch and Trottenbcrg, 1982]. 

Briefly, the multiresolution surface reconstruction algorithm features (i) multiple representations 
of surface shape over a range of spatial resolutions, (ii) local, iterative (relaxation) processes 
that propagate smoothness constraints within each representational level, (iii) local coarse-to-fine 
(prolongation) processes that allow coarser representations to constrain finer ones, (iv) fine-to-coarse 
(restriction) processes that allow finer representations to constrain and improve the accuracy coarser 
ones, and (v) a multilevel coordination strategy that enables the hierarchy of representations and 
component processes to cooperate towards increasing the computational efficiency, usually by orders 
of magnitude. 

Fig. 13 depicts the structure of the algorithm schematically. In this particular case, only three 
levels are shown. Note the 2:1 resolution reduction between adjacent levels. Not only does this 
ratio simplify the component processes considerably, but it is also nearly optimal with regard to 
total computation to convergence (this is conveniently measured in machine independent work 
units, where a work unit is the amount of computation required for a relaxation iteration on the 
finest level) [Brandt, 1977]. The diagram illustrates the intralevel relaxation processes, as well as the 
fine-to-coarse restriction and coarse-to-fine prolongation processes that communicate between levels. 
The figure shows synthetically generated scattered orientation and depth constraints consistent with 
a hemispherical surface. The algorithm reconstructs a dense representation of surface at three 
resolutions. The sparse information at any particular scale can be thought of as a set of constraints 
which defines a discrete surface approximation problem at that level. It is natural then to view the 
multiresolution surface reconstruction algorithm as iteratively solving a coupled hierarchy of discrete 
surface reconstruction problems. 3 For a detailed description of the algorithm see [Terzopoulos, 
1982, 1983a]. 


6. Experimental Analysis of the Algorithm 


The multiresolution visible-surface reconstruction algorithm was tested on a variety of data sets 
including synthetic data, structured light (laser) range data, automated stercopsis and photometric 
stereo data from natural images, and digital terrain model data. Some results are presented 
in this section (for further details and examples, see [Terzopoulos, 1984]). In all the examples 

2 It is possible to accelerate the basic relaxation methods so that fewer iterations are required. However, prac¬ 
tical accelerated methods such as the conjugate gradient method, successive overrelaxation, and Chebyshev 
semi-iteration use global procedures to determine the acceleration parameters. In parallel implementation, 
the greater complexity of the globally accelerated methods and, even more importantly, the communications 
costs of performing the global operations nullifies any potential gains. 

3 A recursive multilevel coordination strategy was employed in the experiments described next. The recursive 
strategy activates only a single level at any one time. We have recently developed a concurrent strategy 
based on a multilevel variational principle [Terzopoulos, 1985b]. Concurrent coordination maintains all 
levels active simultaneously, thus achieving full parallelism. 
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coarse medium fine 




Figure 13. The structure of the multiresolution surface reconstruction algorithm. Iterative relaxation processes 
operate at each level. Fine-to-coarsc and coarse-to-fine processes transfer information between levels. Synthetic 
orientation and depth constraints input to the algorithm are shown at the top. The dense multiscale surface 
representation is output at the bottom. 
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presented, die intralcvcl process was Gauss-Scidel relaxation and die algorithm was started from 
zero initial approximations on all the levels. The border nodes on each grid were preset as depth 
discontinuity nodes to introduce natural boundary conditions which free the reconstructed surface 
on the boundary of fi. 

6.1. Synthetic Data 

The first two examples involve randomly placed depth constraints. The left half of Fig. 14 
shows 15%-density constraints at three resolutions. These constraints were obtained by sampling 
a hemisphere whose 2 values were multiplied by a radial sinusoid. The nodes outside the circular 
region occupied by the constraints were specified as depth discontinuities. The reconstructed surface 
representation is shown on the right half of Fig. 14. In Fig. 15, the 15%-dcnsity depth constraints 
shown on the left are samples of a stacked set of planar surfaces at three resolutions. In this 
example, depth discontinuities were placed along the circular arcs bounding the planes, and along 
the outer edges of the grids. ITie reconstructed surface representation is shown on the right half of 
Fig. 15. This example indicates that discontinuities can be placed along arbitrary contours within fi 
to prevent surface shape from being degraded by unwanted smoothing over sharp depth changes. 

The next examples involve reconstructions from orientation constraints. The left half of Fig. 16 
shows in perspective a set of orientation constraints over a square region. On each of three scales, 
the region is divided into four quadrants each containing constant orientation constraints, and the 
nodes along their boundaries are preset as orientation discontinuities. The surfaces reconstructed 
by the three-level algorithm are shown on the right half of the figure. Since absolute depth cannot 
be determined solely from orientation constraints, a relative depth reconstruction results, with the 
center of gravity of the resulting pyramidal surface resting near the x-y plane. 

The left half of Fig. 17 shows 30%-density scattered orientation constraints consistent with a 
hemispherical surface at three resolutions. The reconstructed surface representation is shown on 
the right. All nodes outside die hemispherical surface patch were specified as dcpdi discontinuities. 
Again, the center of gravity of die surface rests near the x-y plane. 

The next examples demonstrate die integration of both depth and orientation constraints. The 
left half of Fig. 18 shows 15%-dcnsity dcpdi constraints consistent with a hemispherical surface at 
three resolutions. On die right are 15%-density orientation constraints consistent with the same 
surface. Nodes outside die surface have been specified as depth discontinuities. The reconstructed 
surface is shown in Fig. 19. Whereas in the previous example (Fig. 17) only relative depdi can be 
determined for lack of any depth constraints, in the present example the additional depth constraints 
enable the absolute depth of the surface to be determined at all points, hence the surface is “raised” 
to the correct height above the base plane. In addition, note that (10%) uniformly distributed 
noise has been added to the constraint values. Widi the given constraint parameters, die surface is 
slightly bumpy on the finest level. This can be reduced by decreasing the constraint parameters, in 
effect, loosening the springs of the physical model. 

6.2. Structured Light Data 

The multircsolution algorithm was applied to die reconstruction of several objects from raw 
range data supplied by a laser rangefinder constructed by P. Brou at MIT. The scan resolution 
in die y direction is half that in die x direction. A four-level surface reconstruction algorithm 
was employed in die examples. The data was introduced as depth constraints at die finest level 
and transferred to die coarser levels by successive 2x2 averaging between levels. To expediently 
segment the objects from the background, values smaller dian a threshold were treated as depdi 
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Figure 15. Reconstruction of planes with circular depth discontinuities from depth constraints. (Grid 
dimensions: N } f x Njj' = 22 x 17, Nj}* x = 43 x 33, x JV* 3 = 85 x 65. Grid spacings: 

hi = 0.4, h 2 = 0.2, hz = 0.1. Constraint parameters: o<d h] — 2.0 fhj. Computation: 20.375 work units.) 
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Figure 16. Reconstruction of a pyramidal surface with orientation discontinuities from orientation constraints. 
(Grid dimensions: N z l = Njj* - 17, N' z l 2 - N£ 3 = 33, N'j 3 = Njj 3 = 65. Grid spacings: hi = 0.4, 

hi — 0.2, hz = 0.1. Constraint parameters: n p 3 = a q 3 = A. 0 /hj. Computation: 19.5 work units.) 
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Figure 17. Reconstruction of a hemispherical surface from scattered orientation constraints. (Grid dimensions: 
N = N = 17, = Ny ' J = 33, '= C5. Grid spacings: h x = 0.4, h 2 = 0.2, = 0.1. 

Constraint parameters: «p ; = = A.Q/hj. Computation: 22.125 work units.) 
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Figure 18. Depth constraints (left) and orientation constraints (right) consistent with a hemisphere at three 
resolutions. 
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Figure 19. Reconstruction from depth and orientation constraints in Fig. 18. (Grid dimensions: Nl/ 1 = 
Nh' = 17. = 33, N^ = N^ --- 65. Grid spacings: h x = 0.4, h 2 = 0.2, h z = 0.1. 

Constraint parameters n ( i hj = 2.0/hj, a,, j = = 4.0/fcy. Computation: 17.75 work units.) 
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Figure 20. Reconstruction of a lightbuib from range data. (Finest grid dimensions: x Ny* — 257 x 281. 

Grid spacings: hi = 0.8, h 2 = 0.4, /i 3 = 0.2, and /i 4 = 0.1. Constraint parameters: a/j = 0.2/hj. 
Computation: 9.78 work units.) 


discontinuities. Fig. 20 shows the reconstructed surface of a lightbuib. The algorithm smoothes the 
noise in the data and reconstructs the missing points. 

6.3. Natural Image Data 

In this section, we apply the multircsolution surface reconstruction algorithm to depth data orig¬ 
inating from natural images. The examples involve photometric stereo, and two binocular stereo 
algorithms applied to terrain stcrcopairs. 

6.3.1. Photometric Stereo Data 

Photometric stereo is a technique that uses multiple (usually 3) images of a scene from the same 
viewpoint, but with differing illumination [Woodharn, 1981]. Assuming that the surface material 
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Figure 21. Image of a matte white torus. 


is known and that the viewer and lightsources arc far from the object, the method determines 
the surface orientation from the image irradiance. Our surface reconstruction algorithm provides 
a noise resistant technique for computing depth from the surface orientation data provided by 
photometric stereo. We demonstrate this with the image of torus in Fig. 21. The photometric stereo 
data (generated by a system implemented at MIT by K. Ikeuchi) was introduced as orientation 
constraints on a two-level algorithm. Aside from sporadic missing data, the constraints on die 
coarse level arc dense, whereas only every other node on the fine level is a constraint. Fig. 22 
shows the orientation data and the reconstructed torus. 

Our method for reconstructing surfaces from scattered orientation constraints can be compared 
to a variational scheme for obtaining relative depth from dense surface gradient information reported 
by Horn and Brooks [1985j. Their proposed least squares integral f f(v x - p) 2 + (v y - q) 2 dxdy 
will be recognized as being a continuous version of the orientation constraint penalty functional. By 
virtue of the additional smoothness functional S pr {v ), however, our surface reconstruction algorithm 
can deal with orientation constraints that arc scattered. It also can integrate depth constraints from 
other sources to arrive at absolute surface depth. 

6.3.2. Correlation Based Stereo Data 

At the top of Fig. 23 is a stcrcopair on which Kass’s [1983] correlation based stereo algorithm 
was run. The output of the stereo algorithm is shown on the lower left, with brightness proportional 
to disparity. The algorithm has failed to produce a match in the neutral grey patches, so disparity 
is unknown in these areas. To apply the multiresolution algorithm, the disparity data on the finest 
level were reduced by factors of two, through averaging, to three coarser levels. Relatively small 
constraint parameter values were chosen in order to counteract the potentially detrimental effects 
of false matches and noise in the disparity data. The reconstructions on the three coarsest levels 
arc shown as 3D plots in Fig. 24 (the finest level was too dense to represent this way). Fig. 25 
shows isoclcvation contour maps of the solution on all levels. 


1 
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Figure 22. Reconstruction of the torus (right) from the orientation constraints provided by photometric stereo 
(left). (Grid dimensions N = 51 and = 101. Constraint parameters: a h * = 4.0/fy. 

Computation: 52.0 work units.) 


6.3.3. Feature Based Stereo Data 

The next example involves disparity constraints generated by the MPG stereo algorithm 
[Grimson, 1985]. A three-channel version of the stereo algorithm was run on the stcrcopair at 
the top of Fig. 26. The output of the stereo algorithm is shown on the lower part of the figure. 
Disparity information is provided only along zero crossing contours at the three finest scales. In 
the figure, the darkness along contours is proportional to disparity. This disparity data was input to 
a four-level surface reconstruction algorithm. 'The constraints on the coarsest level were derived by 
averaging the constraints from the next finer level. The reconstructions on the three coarsest levels 
arc shown as 3D plots in Fig. 27. Fig. 28 shows isoclcvation contour maps of the solution on all 
levels. 


6.4. Digital Terrain Map Data 
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Figure 23. Natural terrain stereopair (top) and output of Kass’ stereo algorithm (bottom). Ihe images were 
25G x 25G pixels, quantized to 256 levels (provided by the US Defense Mapping Agency). 
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Figure 24. Reconstruction of terrain in Fig. 23. (Grid dimensions: iV'* 1 = JV£‘ = 33, N% 2 = Ny ' J = 65, 
N h 3 = N h 3 = 129t N h< = N h< _ 257 . Grid spacings: - 0.8, /i 2 = 0.4, /i 3 = 0.2, h 4 = 0.1. 

Constraint parameters: - 0.0J//ij. Computation: 29.0 work units.) 
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Figure 25. Isoelevation contour maps of the reconstructed terrain in Fig. 24. 


A four-level surface, reconstruction algorithm was applied to contoured terrain elevation data. 
A contour map of the lllack River Gorges (published by the UK Ministry of Defense) was 
digitized manually on a digitizing tablet by J. Mahoney. The 256 x 256 digital contour array is 
shown at the top of Fig. 29. The constraints input to the algorithm are shown at the bottom 
of the figure. The elevation of the contours is proportional to brightness. Local averaging was 
used to derive the constraints on the coarser grids from those on the finest grid. The terrain 
reconstructions on the three coarsest levels are shown as 3D plots in Fig. 30 (the finest level is too 
dense to represent this way). Fig. 31 shows isoelevation contour plots of the reconstructed terrain 
on all levels. The reconstructed contours on the finest level can be compared subjectively with 
the digitized contours in Fig. 29, but note the reconstructed contours depict elevations half way 
between the original constraint contours for an unbiased comparison. The reconstructed contours 
are somewhat smoother than the (predigitized) contours in the original map — the jaggedness 
introduced by manual digitization has been reduced. The extent of the smoothing can be regulated 
by adjusting the constraint parameters. Shaded image renditions of the reconstructed terrain using 
reflectance map techniques for hill shading [Horn, 1981] are shown at the bottom of Fig. 31. Terrain 
reconstructions using the thin plate surface under tension model were compared to reconstructions 
using the simpler membrane spline model (Laplacian smoothing). The former gives good results, 
whereas the latter generally suffers from insufficient smoothness and produces flat spots across 
terrain peaks [Terzopoulos, 1984] (see also [Bolondi et al., 1976]). 
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Figure 26. Natural terrain stereopair (top) and output of the MPG stereo algorithm (bottom). The images 
were 512 x 512 pixels, quantized to 256 levels (provided by the US Army Engineer Topographic Labs). 
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Figure 27. Reconstruction of data in Fig. 26 on the three coarsest levels. (Grid dimensions: N^ 1 =-- N y l = 33, 
= N y * = 65, N^ 3 = Ny 3 = 129, = Njf = 257. Grid spacings: h L = 0.8, h 2 = 0.4, h z = 0.2, 

hi = 0.1. Constraint parameters: a,i h > = 0.01 /h*. Computation: 31.0 work units.) 
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Figure 28. Isoelevation contour maps of the reconstructed terrain in Fig. 27. 


6.5. Discontinuity Detection Experiments 

The foregoing examples have shown that the surface reconstruction algorithm can handle disconti¬ 
nuities that are prcspecified. In the next two sections, we present examples involving the automatic 
detection of discontinuities. 

6.5.1. The Regularization Approach 

The aerial view stcrcopair of Fig. 32 was input to the MPG stereo algorithm which generated 
the sparse disparity map shown on the top of Fig. 33. The finest level dense disparity map generated 
by a four-level surface reconstruction algorithm is shown at the iower left of the figure. Darkness 
is proportional to disparity. The discontinuities found from this disparity map, using a disparity 
limit G{j > t d = 1 arc shown at the lower right as white contours After the detected points are 
added to the discontinuity map, the surface reconstruction algorithm continues iterating from the 
tentative approximation on the left. The amount of additional computation required is relatively 
small, since the tentative surface is a fairly good approximation in most places. At convergence, the 
reconstructed surface has fractured along the contours to give the solution on the right. Portions 
of the main discontinuities around the buildings have been found, but contours are broken and 
shifted. 

The next example involves the synthesized random dot stereogram in Fig. 7. 'Ihc depth 
constraints generated by a three-channel version of die MPG stereo algorithm are shown in Fig. 34 



















r 
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Figure 29. Digitized contour data (top) and constraints (bottom). The patch to the lower right represents a 
lake. 
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Figure 30. Reconstruction of data in Fig. 29. The terrain reconstruction on the three coarsest levels is 
represented as 3D surface plots. (Grid Dimensions: N% 1 — N* 1 — 33, N'f 2 = Ny' J = 05, N= N ~ 
129, N%* = Ny 4 - 257. Grid spacings: h x = 0.8, -- 0.4, /t 3 = 0.2, h 4 = 0.1. Constraint parameters: 

<*/’ = 0.5/hj.) 
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Figure 31. Isoelevation contour map (lop) of the data in Fig. 30 and shaded representations of the reconstructed 
terrain (bottom). 
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Figure 32. Aerial view of a hospital complex. The stereopair was provided by the UBC Faculty of Forestry. 
Images are 320 x 320 pixels. 


(the finest level dimensions are 320 x 320). The constraints on the coarsest level were obtained by 
averaging those on the next finer level. Fig. 35 shows the smooth disparity maps initially computed 
by a four-level surface reconstruction algorithm. Fig. 36 shows the discontinuities detected from 
these maps with tj = 1. The discontinuities have been superimposed onto the final disparity maps 
in Fig. 37. Better performance is observed in this case due to the simpler surface structure, but the 
contours, while mostly intact, are quite ragged. 

In general, not detecting true discontinuities affects surface shape more adversely over larger 
regions than introducing some spurious ones within a continuous surface. Discontinuity points 
are missed by the thresholding operation, and no adjustment of the global limit can be expected 
to produce perfect results. Note, however, that the surface reconstruction algorithm does not 
break down. Rather, the reconstructed surface degrades as it “leaks” through the gaps. The 
discontinuity detection procedure maty be improved by allowing the disparity limit to vary spatially, 
or by modifying it during multiple passes. On the first pass, surface shape is poorest, so a fairly 
conservative limit should be set to reduce die number of false detections. Conservative limits fail 
to detect many discontinuities, but as more discontinuities are identified, surface shape improves 
and limits can be lowered in subsequent passes to find the less prominent discontinuities. 

6.5.2. The Variational Continuity Control Approach 


A multipass scheme is also employed in the variational continuity control approach to efficiently 
obtain good solutions. An example will be used to explain the strategy. Fig. 38 shows depth 
constraints randomly sampled from a set of sloping planes that form discontinuities along their 
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Figure 33. Discontinuities in the aerial stereogram. Disparity contours generated by stereo algorithm (top), 
full disparity map generated by the surface reconstruction algorithm at the finest level (lower left), and detected 
discontinuities superimposed on the disparity map (lower right). 


si m 
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Figure 34. Depth constraints for the random dot stereogram. 


extremities. A single continuous surface can be reconstructed from these constraints as shown, but 
it smooths over the depth discontinuities and rounds out the orientation discontinuities. 

The algorithm finds both kinds of discontinuities and reconstructs a surface which preserves 
them. When the surface smooths through a depth discontinuity, two spurious regions of high 
curvature border the discontinuity. These spurious regions can easily be mistaken for orientation 
discontinuities. To avoid this unwanted interaction which can substantially slow down the optimiza¬ 
tion process, the algorithm postpones the orientation discontinuity detection phase until all depth 
discontinuities have been found. The surface evolves in several steps over which the parameters 
and /3% in (29) are modified. Each step consists of first flipping the value of the continuity 
control parameter (pj*y or r/' •) from 0 to 1 or conversely, if this lowers the energy (27), and then 
running the reconstruction algorithm to convergence (which always results in equilibrium, since the 
variational principle is convex for fixed pfc ■ and rjv). 

For depth discontinuities. is initially set to a high value that heavily penalizes their insertion, 
then lowered in steps. This strategy of least commitment finds the prominent discontinuities earliest, 
improving the surface as it docs so, and leaves the more subtle ones for later. It results in the 
flipping of relatively few variables in each stage, hence the solution is obtained efficiently. Beginning 
with the continuous surface Fig. 39(a), Fig. 39(b-d) illustrates the steps of die evolving discontinuity 
detection process, during which discontinuities arc determined with increasing accuracy as (3% is 
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Figure 35. Full disparity maps without discontinuities. 


lowered. The energy can be lowered further still if ( 3 % is then increased slightly to eliminate 
spurious discontinuities in Fig. 39(d). Note that since the surfaces have now separated, a very large 
increase would he needed to flip a true discontinuity point (a hysteresis effect). The improved 
surface in Fig. 39(c) results. Next, the orientation discontinuity detection phase is activated and it 
runs in the same way, but modifies ( 3 In this example, the orientation discontinuities arc found 
in only one step. 

Fig. 39(f) shows the final solution. The depth and orientation discontinuities have been made 
explicit and are preserved by the reconstructed surface. Incidentally, the global optimum of the 
variational principle has been found in this example; however, this procedure can generally be 
expected to yield good, though not necessarily optimal approximations. Its main attractions are 
that it is deterministic and efficient. 


7. Discussion and Research Directions 


Several issues concerning the framework for computing visible-surface representations arc discussed 
in this section, and directions for future research arc suggested. The discussion focuses on 
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Figure 36. Detected discontinuities. 


discontinuity detection, choosing constraint parameters, handling rivalries in constraints, grouping 
constraints, invariance properties of the surface reconstruction model, and visible-surface analysis. 
[Terzopoulos, 1984, Ch. 11] contains a more extensive treatment of these and other issues, including 
multiresolution relative depth representations of surfaces, and the possibility of computing visible- 
surface representations “instantaneously” by analog networks. 


7.1. On Discontinuity Detection 

Some recent work in image restoration is of relevance to the problem of piecewise continuous 
surface reconstruction. A piecewise, constant image model employed by Blake [1983] for image 
reconstruction is interesting in that it incorporates “weak constraints” which can be broken at a cost. 
The resulting optimization problem is related to our variational continuity control approach, but 
more restricted. Blake used an adaptive method, which he referred to as “graduated nonconvexity,” 
to obtain good solutions to the nonconvcx problem. It has not been established however whether 
this interesting method applies to the sparse data case as well. 

Gcman and Gcman [1985] used Markov random field models with associated Gibbs distributions 
to restore piecewise constant images corrupted by additive Gaussian noise. The restoration seeks 
a maximum a posteriori estimate of the original image, given die degraded image, and includes 
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Figure 37. Discontinuities and final disparity maps. 


an explicit “line process” that estimates the locations of step edges in intensity. This work was 
restricted to dense image data. The Gemans’ approach was adopted with encouraging results to 
surface reconstruction from sparse depth constraints by Marroquin [1984]. His markov random field 
model, while less restrictive than the Gemans' piecewise constant one, in fact models a membrane 
spline whose smoothness is insufficient for computing visible-surface representations. A line process 
essentially equivalent to the Gemans’ was incorporated to estimate depth discontinuities. The 
numerical solution strategy in both of the above studies was stochastic optimization using the 
Metropolis algorithm and simulated annealing to optimize the nonconvcx functional [Kirkpatrick 
el al ., 1983]. T his strategy can find optimal solutions, but for such large reconstruction problems it 
has been observed to converge notoriously slowly. Based on our experience, we believe that it can 
be accelerated, perhaps enough to make it practical, through die use of multircsolution processing. 

Obviously, the line processes used in the above work as well as our own encoding of disconti¬ 
nuity contour configurations is unplcasingly heuristic and in need of refinement. The discontinuity 
map can be augmented by nodal variables to encode the local orientations of the curvilinear ele¬ 
ments to a higher degree of accuracy. Such an encoding is employed by Zuckcr and Parent [1984] 
in an optimization (relaxation labeling) approach to finding contours in images. It appears that 
ideas from their work can also be applied to finding surface discontinuities within our framework. 

A promising possibility is to employ ID controiled-continuity stabilizers as formal models 
of smoothness constraints along surface discontinuity contours in the x-y plane. A functional 
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Figure 38. Scattered depth constraints consistent with sloping planes meeting discontinuous^ (top) and the 
smooth reconstructed surface (bottom). 


that naturally comes to mind is the curvilinear analog of the thin plate surface under tension: 

j3 b (s){ft a (s)(d 2 c/ds 2 ) + [l-j3 a (s)](dc/ds)} ds, where s denotes arc length along discontinuity 
contours c e C. Here, /3 b allows breaks, while /3 a allows angles (tanjcnt discontinuities) to form 
in the discontinuity contours. Again, additional energy penalties must be associated with these 
occurrences. Given our finite element representation of surfaces, curvilinear finite elements are 
the natural local representation for discontinuity contours. The combined variational principle has 
both a surface component and an analogous contour component. Although technically nontrivial, 
a formulation of surface reconstruction generalized along these lines has very strong appeal. 

7.2. On Constraints — Parameters, Rivalries, and Grouping 

The constraint (spring) parameters offer the flexibility to individually tune the coerciveness of each 
constraint on the reconstructed surface. In the special case of Gaussian error distributions, the 
parameters should be inversely proportional to the expected variances (a t = 1/A<r?). It ought to be 
possible for the low-level visual processes to associate a variance estimate or confidence with each 
constraint that they provide. In general, however, it’s not obvious how to choose the constraint 
parameters optimally. 

The constant of proportionality A -1 can also be used to tune the overall smoothness of 
the reconstructed surface. Cross validation techniques may be used to set A optimally (e.g., 
[Wahba and Wendclbcrger, 1980]). The basic criterion is to choose A so as to minimize over all 
constraints the (weighted) discrepancy between each constraint and its value as estimated from the 
surface reconstructed using the remaining constraints. Unfortunately, this involves computationally 
expensive sequential algorithms. Interestingly, the continuous tuning of surface smoothness is 
analogous to the scale space filtering technique proposed by Witkin [1983] with the added attraction 
that it can be applied to scattered data. 

Although the variational principle was designed to account for measurement errors in the 
constraints, the possibility of massive rivalries between constraints from different sources, such as 
stereopsis and analysis of motion processes, was disregarded. Massive rivalries are unnatural visual 
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Figure 39. Evolution of the discontinuity detection process. 


phenomena that can nevertheless occur, especially under contrived conditions, and they often lead 
to multistable percepts [Attneavc, 1971]. The framework can potentially accommodate rivalries with 
a mechanism that inhibits individual or entire sets of constraints by nullifying selected constraint 
parameters. This mechanism can be activated by a global arbitrator which monitors the contents 
of visible-surface representations to detect rivalries may also have access to higher level knowledge 
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about the scene. "Die arbitrator’s influence can account for multistability. 

A particular type of rivalry arises from transparent surfaces. For instance, a surface such as a 
dirty window in front of a background scene would lead to two well defined populations of depth 
(and orientation) constraints over the same visual angle, one from the window, the other from the 
background. A transparency interpretation can be arrived at by an arbitrator which monitors the 
surface reconstruction process looking for high approximation error between surface and constraints 
over a significant area. Under these conditions the arbitrator can trigger a constraint grouping 
process which clusters the constraints into two populations, based on depth values, say. Multiple 
surfaces can then be reconstructed over the same visual area for each resulting constraint population. 
This scheme has been applied on a transparent surface random dot stereogram [Terzopoulos, 1984, 
Ch. 11]. 

7.3. On Invariance Properties of the Surface Model 


As a transformation from sparse constraints to dense surfaces, the thin plate under tension model 
can be shown to be invariant under (i.e. Commutes with) certain image plane transformations 
applied to the constraints; namely, translations, rotations, and similarity transformations. This 
implies that surface shapes will be preserved through rigid motions of the scene or viewpoint 
parallel to the image plane or along the view direction. These are essential invariance properties 
for visible-surface reconstruction [Terzopoulos, 1982]. 

Note, however, that the thin plate spline, characterized by the small deflection approximation 
//«4 + 2 vl y -I- v yy dx dy to the bending energy density of a thin plate, is not invariant under 
arbitrary 3D transformations of the constraints. Thus, surface interpolation using this expression is 
not invariant under changes in the view direction, as Blake [1984] points out. He shows that rotating 
the view direction induces the ID analog of the thin plate spline to “wobble,” and he demonstrates 
that this effect is most pronounced as the (continuous) spline is inclined sharply with respect to the 
viewer or is forced to bend sharply. Blake views this as a problem that should be eliminated by 
employing the large deflection bending energy of the thin plate, a convex combination of the mean 
and Gaussian curvatures of the surface v(x,y) f which is view direction invariant. 

Although £ pr {v) can also be made view direction invariant by employing the large deflection 
counterparts for the thin plate and membrane bending energies, this approach has a serious 
technical drawback [Terzopoulos, 1984]: The large deflection formulas lead to an extremely difficult 
nonlinear problem (e.g., the large deflection equations for the thin plate are two coupled nonlinear 
fourth-order partial differential equations known as Von Karmann’s equations [Szilard, 1974]). 

Fortunately, the surface reconstruction model, as it stands, is not hampered by die lack of 
view direction invariance because the available constraints are usually sufficiently dense in practice 
to tightly determine surface shape; as the view direction is varied, the reconstructed surface would 
vary negligibly (note that Blake’s experiments reveal a significant wobble effect just in die case 
of extremely sparse constraints). Furdiennore, the explicit introduction of depth and orientation 
discontinuities alleviates much of the wobble precisely at those places where Blake’s experiments 
show it to be most pronounced on a globally continuous surface. An interesting psychophysical 
experiment would be to determine whether there might be some slight variance in the surfaces 
perceived by humans viewing sparse random dot stereograms while the dots undergo simulated 
rigid 3D transformations and, if so, whcdier die variations arc consistent with die reconstruction 
model (J. Mayhcw, personal communication). 
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7.4. On Visible-Surface Analysis 

The visible-surface representation is an intermediate and volatile description of the 3D surfaces 
in scenes. It drives ensuing processes which generate stable higher-level representations of shape 
that are better tuned to object recognition. The processing begins with visible-surface analysis 
whose goal is to abstract from the numeric, viewer-centered representation a rich set of more 
symbolic, object-centered features that are stable through viewpoint changes. The extraction of 
geometric surface features is facilitated by the dense shape information provided by visible-surface 
representations. 

A promising approach to visible-surface analysis is to apply concepts from differential geometry 
[do Carmo, 1976]. For instance, a surface’s intrinsic geometry (including Gaussian curvature, 
geodesics, etc.) is determined completely by the first fundamental form, which defines arc length 
over the surface. It’s extrinsic geometry (including normal curvature, principal curvatures, etc.) are 
determined by the second fundamental form, which describes the deviation of the surface from the 
local tangent plane. The fundamental theorem of the local theory of surfaces (usually attributed 
to Bonnet) states that the analytic study of surface properties consists of the study of the two 
fundamental forms; i.e., the six fundamental tensor coefficients (which arc not all independent) as 
functions of the two independent parameters of the surface. The fundamental forms are invariant 
under changes in the parameterization, and together they determine surface shape up to rigid body 
transformations. These properties make them ideal foundations for object centered symbolic surface 
representations. 

The visible-surface representation makes it possible to estimate the first and second funda¬ 
mental forms on a point-by-point basis over the entire visible surface. The finite element shape 
representation reduces the computation of crucial local surface features such as the Gaussian curva¬ 
ture, principal curvatures, and principal directions to the evaluation of simple algebraic expressions 
of neighboring nodal variables (see [Terzopoulos, 1984, Ch. 11] for derivations), ft is then a simple 
step to determine the elliptic, hyperbolic, parabolic, umbilic, and planar points, as well as geodesics, 
asymptotes, and lines of curvature. 

For example, Fig. 39 shows the reconstructed surface of a lightbulb. Fig. 40 shows the Gaussian 
curvature K ( x , y) computed for the reconstructed lightbulb surface of Fig. 20. The elliptic points 
(K > 0) are shown in white, the hyperbolic points (K < 0) are shown in black, and the parabolic 
(K = 0) points separate the two regions. Note the alternation in the sign of curvature at the 
screw mount. Fig. 41 plots the computed field of principal directions for the lightbulb at the 
two coarsest scales. These demonstrations illustrate the feasibility of reliably computing from these 
representations higher-order intrinsic and extrinsic properties of surface shape. The reliability can 
be attributed to the regularizing properties of the thin plate surface under tension which overcomes 
the potentially detrimental effects of noise in the data, while preserving discontinuities. For further 
analysis of the kinds of features that can be computed from dense, numeric, representations of 
surfaces see, e.g., [Brady et al ., 1985] or [Medioni and Nevada, 1984]. 


8. Conclusion 


Constraints on surface shape, contributed by multiple low-level visual processes, can be computed 
reliably at multiple resolutions, but only at scattered locations in the field of view. Subsequent 
visual processing can be facilitated substantially if the scattered constraints are transformed into 
visible-surface representations that make surface shape explicit everywhere. To accomplish this 
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Figure 40. Elliptic (white) and hyperbolic (black) points of the reconstructed lightbulb at four scales. 


effectively, information must be integrated over multiple visual modalities and fused across multiple 
scales of resolution. 

In this paper, we have developed a computational theory of visible-surface representations. 
Within a unified computational framework, formal solutions were offered to fundamental problems 
of reconstructing visible surfaces: (i) integrating constraints on the depth and orientation of surfaces 
across various modalities and scales, (ii) interpolating surface shape information into (piecewise) 
smooth surfaces, (iii) discovering discontinuities in surface depth and orientation and enabling them 
to restrict interpolation, and (iv) efficiently maintaining consistency in distributed, multi resolution 
visible-surface representations. 

A visible-surface reconstruction algorithm implements the framework. Extensive testing has 
shown it to be viable. The algorithm coordinates cooperative processes within a multircsolution 
hierarchy of surface representations to dramatically increase computational efficiency. It is well 
suited to implementation on massively parallel networks of simple, locally interconnected processors. 
Such computational networks are suggestive of biological mechanisms and are also well suited to 
VLSI technology. 
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Figure 41. Principal directions for the reconstructed lightbulb at the two coarsest scales. The directions of 
greatest curvature are shown on the left. Those of least curvature are shown on the right. 
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